Embedding approach for dynamical mean field theory of strongly correlated 

heterostructures 
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We present an embedding approach based on localized basis functions which permits an efficient 
application of the dynamical mean field theory (DMFT) to inhomogeneous correlated materials, such 
as semi-infinite surfaces and heterostructures. In this scheme, the semi-infinite substrate leads con- 
nected to both sides of the central region of interest are represented via complex, energy-dependent 
embedding potentials that incorporate one-electron as well as many-body effects within the sub- 
strates. As a result, the number of layers which must be treated explicitly in the layer-coupled 
DMFT equation is greatly reduced. To illustrate the usefulness of this approach, we present nu- 
merical results for strongly correlated surfaces, interfaces, and heterostructures of the single-band 
Hubbard model. 

PACS numbers: 73.20.-r, 71.27. +a, 71.10.Fd 



I. INTRODUCTION 

In recent years, there is growing interest in the elec- 
tronic properties of surfaces and interfaces of strongly 
correlated materials^ For instance, the discrepancy be- 
tween photoemission spectra and theoretically derived 
bulk spectra of a number of transition-metal oxides has 
been attributed to changes in the electronic structure at 
the surface of these materials Regarding the inter- 
face, heterostructures made out of thin atomic layers of 
perovskite-type oxides have been the target of intense 
study as promising candidates for electron-correlation- 
based devices ji^i^ A well-known example is the inter- 
face between LaTi03 and SrTiC>3, which exhibits metal- 
lic behavior in spite of the fact that the two constituent 
bulk materials are insulators, 5 

On the theoretical side, inhomogeneous layered sys- 
tems have been studied by several authors within the dy- 
namical mean field theory^ (DMFT), which contributed 
significantly in the last decade to the understanding of 
a variety of strongly correlated bulk materials J^- Pott- 
hoff and Nolting investigated the metal-insulator tran- 
sition (MIT) at the surface of the single-band Hub- 
bard modeL 12 ' 13 Liebsch studied the t% g valence bands of 
perovskite-type oxides such as SrV03 by using a three- 
band tight-binding Hamiltonian and showed that elec- 
trons at the surface are more strongly correlated than 
in the bulk due to the reduction of the effective sur- 
face band widthJ^ Helmes et at. considered a metal- 
insulator interface within the single-band Hubbard model 
and studied the scaling behavior of the metallic pen- 
etration depth into the Mott insulator near the criti- 
cal Coulomb energy^ Okamoto and Millis investigated 
the electronic structure of heterostructures in which a 
finite number of Mott-insulator layers were sandwiched 
between band insulators i 16 ' 17 Analogous calculations for 
heterostructures consisting of correlated model systems 
were also carried out by Kancharla and Dagottoi& and 
Riiegg et at 19 Modulation doping effects at heterojunc- 



tions were investigated by Oka and Nagaosa3£, Lee and 
MacDonald^i and Gonzalez et alM Electron transport 
through a nano-size correlated-electron system connected 
to metal electrodes was studied by combining DMFT and 
a non-equilibrium Green-function techniquei 23 ' 24 

To solve the DMFT equation for inhomogeneous lay- 
ered systems one needs to construct the lattice Green 
function of surfaces or interfaces consisting of an infinite 
number of atomic layers. While this is feasible within a 
linearized version of DMFT, 12 for a complete numerical 
solution of the DMFT equation most previous calcula- 
tions employed a slab model consisting of a finite num- 
ber of layers to simulate the system. Although finite- 
size effects can be reduced by systematically increasing 
the number of layers, the one-electron density of states 
(DOS) projected on each layer converges rather slowly 
with increasing number of layers since the energy levels 
in the normal direction are discrete. Hence, it is desir- 
able to develop a method for solving the DMFT equa- 
tion for truly semi-infinite surfaces and interfaces be- 
tween two semi-infinite materials. Chen and Freericks 
solved the DMFT equation for a thin doped Mott in- 
sulator sandwiched between two semi-infinite metals by 
applying the quantum zipper algorithm to the Falikov- 
Kimball Hamiltonian. 25 In the present work we pursue 
a different approach by extending the concept of tight- 
binding embedding, originally developed for the evalua- 
tion of the electronic properties of defects in solids, to 
the DMFT for inhomogeneous layered systems. 

We employ a localized basis set to describe the Hamil- 
tonian of the system. The heterostructure is divided 
into a central interface region containing a finite num- 
ber of atomic layers, f2, and two adjacent semi-infinite 
bulk regions coupled to Q. The interface region is as- 
sumed to include the first few surface layers of the ac- 
tual substrates. Both the central region and the sub- 
strates may exhibit strong correlation effects. Within 
the one-electron approximation, the effects of an adja- 
cent semi-infinite system on f2 can be expressed as a 
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complex energy-dependent potential acting on the Hamil- 
tonian matrix of O, which is called "tight- binding em- 
bedding potential" The same quantity is called 
"contact self-energy" in transport theory based on the 
non-equilibrium Green function formalisrm 29 i 30 Here, we 
extend this embedding approach in order to include 
Coulomb correlations in the substrate within the single- 
site DMFT. Thus, the energy-dependent embedding po- 
tential accounts for one-electron and many-body effects 
within the substrates. The advantage of this extension is 
that the layer-coupled DMFT equation for a non-periodic 
surface or interface system made up of an infinite number 
of atomic layers is greatly simplified since only a small 
number of layers belonging to Jl needs to be treated ex- 
plicitly. The embedding potential is derived from a sep- 
arate DMFT calculation for the adjacent bulk systems. 

The plan of this paper is as follows: In Sec. |TT] we 
present our formalism that combines the tight-binding 
embedding theory with DMFT. As examples, in Sec. IIIII 
we present numerical results for surfaces and interfaces 
of the single-band Hubbard model at half filling. A sum- 
mary is given in Sec. IIVI 

II. THEORY 

A. Hamiltonian 

We take the z axis as the surface normal pointing from 
left to right. The j th atomic layer is located at z = zj 
(zj > Zj-i). The position of the m th atom in layer j 
is denoted by r a — (x„ y ,Zj), where the index a repre- 
sents a pair of indices (m,j). The localized basis function 
centered at r a with orbital index (3 and spin index a is 
denoted by ipa(3cr- The basis set {i/japo-} is assumed to be 
orthonormal. Hereafter, we use indices with tilde such as 
p and q to refer to basis functions in the basis set {ipa/3a}- 
With this abbreviated notation, the one-electron part of 
the Hamiltonian is written as 

h = *m 4 C «> w 

where ct and Cq are the creation and annihilation oper- 
ators, respectively, and summation is taken over pairs 
(p, q) having the same spin and located on the same 
or nearby atomic sites. The one-electron Hamiltonian 
may be derived, for example, from a first-principles 
electronic-structure calculation within density-functional 
theory through the use of maximally localized Wannier 
functions M- 

In the present work we consider onsite Coulomb inter- 
actions, 

V = 2 ^ ] Vpq,fs C p C q C s c f ? (2) 
p,q,r,s 

where p, q, f, and s are located on the same site, and in 



addition, p and f (q and s) have the same spin. The full 
Hamitonian of the system is given by H = h + v. 

We now divide the system into three parts. The cen- 
tral region with atomic layer index j running from 1 to N 
is called Q. The semi-infinite region with layer number 
j < 1 is called "left substrate" L, and the semi-infinite 
region with j > N is called "right substrate" R. In the 
case of a semi-infinite surface, it is understood that the 
system consists only of £1 and R. In the following, we 
present the theory for the interface geometry. The anal- 
ogous equations theory for a semi-infinite surface are de- 
rived straightforwardly by omitting all terms with index 
L. The one-electron Hamiltonian h in Eq. ([T|) is decom- 
posed into seven parts: 

h = fiLL + hnn + Hrr + ^ (j>>QJ + hjo) , (3) 

J=LM 

with 

h JK = J2(J\h\k)c]c- k , (4) 

where J and K denote one of the three regions, L, f2, and 
R, and the basis function j (k) belongs to region J (K). 
It is to be noted that the matrix elements of the inter- 
regional terms, t^j, where w £ fl and j e J ( J = L, i?), 
are non-vanishing only when w and j are close to the 
boundary between O and J, since transfer integrals are 
short-ranged. For the same reason, hm and h R L are 
assumed to vanish. 



B. Non-interacting Green function 

As a brief review of the tight-binding embedding 
theory ) 26 i 27 i 28 we outline first the calculation of the Green 
function (resolvent) of the one-electron Hamiltonian h, 

gp q (e) = (p\(e-h)- 1 \q). (5) 

When both indices of this Green function belong to SI, the 
tight-binding embedding theory reveals that the interac- 
tion with the left and right substrates can be expressed 
in terms of embedding potentials acting on Jl: 

4*'(<0 = f*r4( e ) W> (6) 

where w,w' G S7 and the summation is implied for re- 
peated indices. g L and g R are the Green functions of the 
left and right substrates, respectively, and are defined as 

9 L (e) = {e-h LL )-\ (8) 
g R (e) = (e-hnn)- 1 . (9) 

It should be noted that s L and s R in Eqs. © and (JT]) 
are non-vanishing only when both w and w' are located 
close to the boundaries of fL 
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Using these embedding potentials, the Green function 
defined in Eq. ([5]), when both indices belong to f2, can 
be calculated as 



Sto'(e) = (w\[e- h em (e)] 1 | w'), 



(10) 



where the effective Hamiltonian in the embedded region 
n is given by 



L = hnn + s L (e) + s R (e)- 



(11) 



Thus, the calculation of a system which is non-periodic 
in the z direction is reduced to the inversion of a matrix 
defined in O with a finite thickness. 

We point out that, inspite of the finite size of the cen- 
tral interface region, the use of the complex embedding 
potentials ensures that the spectral distribution is con- 
tinuous. In particular, for a uniform system with layer- 
independent energy levels and hopping matrix elements, 
the local DOS of each layer coincides with the bulk DOS. 
Thus, there are no discretization effects stemming from 
the finite number of layers in the central region. 



C. Dynamical mean field theory 

We now incorporate the Coulomb interactions and cal- 
culate the finite-temperature Green function of the full 
Hamiltonian, H — h + v. The effects of the Coulomb 
interactions can be described by a frequency-dependent 
self-energy E(iw„), where u n = (2n + \)ttT are Mat- 
subara frequencies at temperature T. In the present 
work, we restrict ourselves to the single-site approxima- 
tion and assume that the matrix elements Y>pq(ioj n ) are 
non- vanishing only when p and q are on the same atomic 
site. Hence, 



£(iu>n) = X] ^JK(i^n)SjK, 



(12) 



J,K 



where Tijk is defined in the same way as Eq. (|4|) with h 
replaced by E. The lattice Green function of the whole 
system is defined by 

Gpg(iuj n ) = (p\[iu n + pb-h- 'E(iuj n )]~ 1 \q), (13) 

where /i denotes the chemical potential of the system. 

As in the case of non-interacting systems, we define 
the embedding potentials of the left and right correlated 
substrates as 



where G L {iuj n ) and G R (iu n ) are defined by 



G L (iu; n ) 
G R (iuj n ) 



A* - h LL 
A* - h RR 



' ^RR^n) 



(14) 
(15) 



(16) 
(17) 



With these definitions, the lattice Green function defined 
by Eq. (|13|). when both indices belong to f2, is calculated 



Gww>(iu n ) = {w\[iuj„ + At - H cm (iu n )] \w), (18) 
where the interacting embedded Hamiltonian is given by 
H em = hnn + Snn(iw n ) + S L (iu n ) + S R (iu>„). (19) 

Suppose now that both substrates are semi-infinite 
crystals having three-dimensional translational symme- 
try and that the boundary between f2 and L (i?) is posi- 
tioned a few atomic layers toward the interior of the crys- 
tal such that the electronic structure in L (R) converges 
to that of the bulk. We may then assume that the matrix 
elements of Ejj on all atomic sites in J (J = L, R) be- 
come identical with those of the Coulomb self-energy on 
the corresponding atomic site in the bulk crystal. There- 
fore, we are left with determining the self-energy in the 
embedded region, Sjiq- This can be achieved via the 
following three steps: (i) perform a standard DMFT cal- 
culation for the bulk crystals corresponding to the left 
and right substrates to obtain the Coulomb self-energies 
in the bulk, (ii) construct the embedding potentials of 
both substrates, S L and S R , and (hi) perform a layer- 
coupled DMFT calculation in the embedded region Q to 
self-consistently determine Snn- 

The embedded DMFT calculation in the third step is 
conducted in a standard manner. Starting with an input 
lattice self-energy Enn, one calculates the lattice Green 
function G(iu) n ) in £1 by using Eq. (fTS]) . To avoid double 
counting of local Coulomb interactions, it is necessary to 
remove at each atomic site in J7, r a , the onsite Coulomb 
self-energy term from the lattice Green function. This 
yields the bath Green function, 



G ,aa{iuJn) = G aa (iuj n ) 1 + Y, aa (iuj n ) 



(20) 



where E aa is the projection of Eqq on atomic site r a 
defined by 

E QQ (iw„) = ^ (ui 1 Soo(i^n) | w'^csj', (21) 



with both w and w' located on r a . G aa , the projection of 
the lattice Green function G on atomic site r a , is defined 
in the same way. Both E Qa and G aa are iV a x N a matrices, 
where N a is the number of basis functions centered at r a . 
Within the single-site approximation, Enn is diagonal 
with respect to atomic sites, so that 



(22) 



The bath Green function Go,aa determines the Weiss 
mean- field Hamiltonian at site r a . One then adds the lo- 
cal Coulomb interactions of the form Eq. © and solves 
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the single-site many-body impurity problem at site r a by 
numerical methods, such as the quantum Monte Carlo 
approach , 32 ! 33 exact-diagonalization 3 -^ (ED), or the nu- 
merical renormalization group method^ The resultant 
impurity Green function, G™ p , is used to derive the out- 
put impurity self-energy via 



(23) 



The key assumption in DMFT is now that this impu- 
rity self-energy is a physically reasonable representation 
of the lattice self-energy. Thus, S aQ («o; rl ) « E"° p (iw n ). 
This self-energy is therefore used as input in Eqs. (fT9|) . 
(f20|) in the next iteration. This procedure is repeated 
until the difference between the input and output self- 
energies becomes sufficiently small for all atomic sites in 
the embedded region f2. 



III. RESULTS AND DISCUSSION 

A. Hubbard model 

To demonstrate the DMFT embedding approach we 
present results for the single-band Hubbard model, 

H = e a h aa - 2J (tab c\ a c ba + H.c.) 

aa {ab)o 
+ J^Uahafhal, (24) 



where n n 



and the summation in the second 



term is taken over nearest neighbor sites. We consider a 
simple cubic lattice with its three principal axes oriented 
along the x, y, and z directions. The interface points in 
the (001) direction. In each layer, all sites are assumed to 
be equivalent (lxl structure). We label the site energy 
of layer j as tj , the Coulomb energy of layer j as Uj , the 
x and y components of the inplane transfer integrals in 
layer j as and t y , and the transfer integral between 
two nearest- neighbor layers, j and k, as t z - k . The Hamil- 
tonian parameters in L (R) represents a particular bulk 
crystal with a single atom in the unit cell. The Hamil- 
tonian parameters in Q approach those of the left-hand 
side (right-hand side) crystal near the boundary to L (R) , 
while they are allowed to deviate from these bulk param- 
eters in the interior of region f2. In the present work, we 
consider only paramagnetic solutions and omit the spin 
index a in the discussion below. 

As an input, one needs the embedding potentials of 
both substrates. Let us consider the left substrate L, 
whose site energy, Coulomb energy, and transfer inte- 
grals are given by cl, Ul, t x L , t y L , and t z L . First, we 
ignore the Coulomb interactions and derive the embed- 
ding potential for non-interacting electrons as defined by 
Eq. Because of translational symmetry in the plane, 
the embedding potential is diagonal with respect to the 
two-dimensional wave vector k = (k x ,k y ) and can be 



expressed as sfj (k, e) by introducing a mixed representa- 
tion with k and layer indices in O, % and j. Here, the 
wave vector is measured in units of the inverse of lattice 
constant, i.e., —n < k x ,k y < tt. For the present nearest- 
neighbor transfer model, the only non- vanishing element 
is (k, e) , which is given as 



J f 1 (k,e) = (tg 1 ) 2 g^(e-E L (k)) 



(25) 



where 



E L (k) = e L - 2t x L cos{k x ) - 2t y L cos{k y ), (26) 



and gjj, with j, j' < 1 denotes the Green function of 
a semi-infinite tight-binding chain with nearest-neighbor 



transfer integral, 
Soven^ 



-t? . According to Kalkstein and 



3ooH = 



(27) 



The embedding potential of region L in the presence 
of Coulomb interactions is obtained by incorporating the 
effects of electron correlations in the bulk crystal into Eq. 
(1251) as 



Sf 1 (k,i« n ) = (*S 1 ) 2 9w(w\ 



(28) 



with 



w = iuj n + fj, - E L (k) - (a\E LL (iu} n )\a}, (29) 



where a is any site in L. On the right-hand side of Eq. 
(|29|) . the Coulomb self-energy in L is determined by a 
bulk DMFT calculation. The embedding potential of the 
right substrate R can be constructed in the same way. 

Using Eq. (fT5)l . the lattice Green function in region £1 
is now calculated as 



G a b(iLOn) 



dk 

X (j|[tw n + H- HenfciWn)]- 1 ]?), (30) 



where a — (m,j), b = (m',j'), 1 < j,j' < N. In 
the mixed representation the embedded Hamiltonian, 
H em (k,iu> n ), is an N x iV matrix, 



(j|#em|/} 



E 3 (k)6 33 , 
Sjj,(k,iuj n ) 



3 33' ' 



S&(k,iw„). 



Here. 



Ejik) = ej - 2t| cos(k x ) - 2t y 3 - cos(ky). 



(31) 



(32) 



The Coulomb self-energy of layer j, Y,j(iu} n ) — 
(a|San(w n )|a) with a = (m, j), is diagonal with respect 
to the layer index and has no dependency on k within the 
single-site approximation. As argued above, in a nearest- 
neighbor tight-binding system, only the embedding po- 
tentials Sii(k,iuj n ) and S§ N (k,i/jj n ) are finite. The 
layer-dependent onsite Green function for a = b = (m, j) 
will be denoted as Gj(iu> n ) and the corresponding bath 
Green function as Goj(iu! n ). 
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B. Numerical results 

We consider first the surface of a semi-infinite Hub- 
bard model having uniform Hamiltonian parameters, i.e., 
Uj = U, ej = -U/2, and tj = t v - = t z jk = t for all lay- 
ers including the surface plane. By choosing the chemi- 
cal potential /i as zero, all layers become half-filled due 
to electron-hole symmetry. For zero temperature, the 
same system was studied by Potthoff and Nolting^ who 
showed that there is a uniform critical Coulomb energy 
U c , at which both bulk and surface simultaneously un- 
dergo a metal-insulator transition. For a complete nu- 
merical solution of the DMFT equation, they adopted a 
slab geometry consisting of 10 to 20 atomic layers rather 
than treating semi-infinite surfaces. 

As impurity solver we employ the finite-temperature 
ED method. Thus, for each layer j, the bath Green func- 
tion Eq. (|20p is projected on to a small cluster consisting 
of a single impurity surrounded by several bath levels. 
Eq. ([20]) is therefore approximated as 

iun + n-Ej-Y . 1 J "*' , (33) 

where Ej represents an impurity level for layer j, Ej^ 
the corresponding bath levels, and Vjk specifies the hy- 
bridization matrix. We use rib — 7 bath orbitals in the 
numerical results presented below. The inclusion of the 
ficticious impurity level Ej provides a more accurate pro- 
jection of Go,j(iu> n ) than for a cluster consisting only of 
bath orbitals. The interacting Green function of the clus- 
ter with onsite Coulomb energy U at finite temperature 
T is derived by calculating the low eigenvalues of the clus- 
ter via the Arnoldi algorithm and applying the Lanczos 
procedure for computing the excited state Green func- 
tion. More details of the ED method can be found in 
Ref. M 

In Fig. Q] we show the calculated quasiparticle weight 
Z of the semi-infinite Hubbard model at T = 0.02 in the 
metallic range U < U c as a function of layer index j, 
where the outermost layer corresponds to j = 1. The 
hopping integral is taken to be t = 1 and defines the en- 
ergy scale. The crosses on the vertical axis indicate the Z 
values of the bulk metal determined by a separate bulk 
DMFT calculation. The embedding potential acts on 
layer N on the right edge of the surface region. Solid dots 
and open circles provide the results obtained for N = 5 
and N = 3 embedded layers, respectively. The excellent 
agreement between the two sets of calculations demon- 
strates that the embedding potential represents correctly 
the one-electron as well as many-body properties of the 
semi-infinite substrate. Moreover, it is evident that one 
needs only a few embedded layers to simulate the semi- 
infinite system. Although the latter point is not crucial 
for single-band model systems, for the calculation of re- 
alistic multi-orbital materials the embedding treatment 
yields a substantial reduction of computer time compared 
with slab calculations in which at least ~10 layers must 




Layer j 



FIG. 1: (color online) Quasiparticle weight Z of semi-infinite 
Hubbard model for simple cubic lattice in the (001) orienta- 
tion as a function of layer index j. Temperature T = 0.02. 
Solid dots and open circles are results with N = 5 and N = 3 
embedded layers, respectively. Crosses on the vertical axis 
indicate the bulk Z corresponding to four values of U. Lines 
are drawn as a guide to the eye. 



be explicitly taken into consideration. 

As can be seen in Fig. 1, the calculated Z exhibits 
an oscillatory behavior near the surface which follows 
from the Friedel oscillations of the layer-dependent den- 
sity of states. In the first layer, Z is smaller than the 
bulk value, implying that electrons at the surface are 
more strongly correlated than in the bulk. As discussed 
by Potthoff and Noltingi2 and Liebschjii this is essen- 
tially a one-electron effect arising from the layer depen- 
dence of the one-electron DOS of the cubic tight-binding 
Hamiltoniani2£ Because of the loss of nearest-neighbor 
sites, the effective band width in the first layer is re- 
duced, so that Coulomb correlations at the surface are 
enhanced. 

Next, we study the interface between two semi-infinite 
Hubbard models. We consider a uniform system with re- 
gard to the transfer integrals, i.e., fj = t v - = tj k = t = 1. 
In the left half-space, we choose the Coulomb energy 
as Ul = 6 to represent a good metal with a relatively 
large quasiparticle weight, while we take a variable, larger 
Coulomb energy Ur in the half-space on the right. Fur- 
thermore, by choosing e 3 as — E/l/2 and —Ur/2 in the 
left and right half-spaces, respectively, and by setting the 
chemical potential as fi = 0, all layers are half-filled. The 
same model was recently investigated by Helmes et al^ 
who used the NRG method as impurity solver. These 
authors focused on the critical range for Ur ~ U c and 
discussed the scaling behavior of the metallic penetra- 
tion depth into the Mott insulator. To reduce finite-size 
effects, a relatively thick slab consisting of ~60 layers was 
used to simulate the interface. Also, to avoid numerical 
difficulties stemming from the energy discretization the 
van Hove singularity of the two-dimensional layer DOS 
was cut off at a finite value. 
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FIG. 2: (color online) Quasiparticle weight Z at the interface 
between two semi-infinite Hubbard models for simple cubic 
lattice in (001) orientation as a function of layer index j. The 
boundary between the two metals is marked by a dashed line. 
Temperature T = 0.02. Solid dots and open circles show 
results with N = 10 and N = 6 embedded layers, respectively. 
Crosses on the vertical axis indicate bulk Z values. Lines are 
drawn as a guide to the eye. 



Figure [2] shows the calculated quasiparticle weight Z 
in the metallic range Ur < U c at T = 0.02 as a func- 
tion of layer index j, which is measured here relative to 
the boundary layer of the left-hand side metal. To de- 
scribe the deviation of the electronic structure from that 
in bulk metal, we incorporate in the embedded region 
N = 10 layers (solid dots), of which the left (right) five 
layers possess Coulomb energy Ul {Ur). For comparison, 
we also show for Ur = 12 (open circles) the result ob- 
tained with only N — 6 embedded layers. The excellent 
agreement between the two sets of calculations corrob- 
orates again the efficiency of the embedding method to 
treat semi-infinite substrates. 

The quasiparticle weight Z of the surface layer of the 
good metal on the left-hand side is seen to be reduced 
whereas at the surface of the poor metal on the right- 
hand side of the boundary plane it is enhanced. Evi- 
dently, the good or bad metallic character of one metal 
spills over into the neighboring metal. In contrast to the 
case of the semi-infinite surface discussed above, this is a 
genuine many-body effect, since the one-electron DOS is 
layer-independent if the Coulomb interaction is switched 
off. The deviation of Z from the bulk value in poor metal 
on the right decreases with the distance from the bound- 
ary plane, which is in accord with the work of Helmes 
et alr^ On the left of the boundary plane, Z is seen to 
be weakly modified with respect to the bulk value only 
in the first two layers. Thus, as a result of better elec- 
tronic screening in the good metal, Z approaches the bulk 
value more rapidly than in the poor metal. Essentially, 
one needs to incorporate only one or two layers in the 
embedded region to describe the interface properties of 
the good metal on the left-hand side. 



FIG. 3: (color online) Quasiparticle weight Z of the het- 
erostructure consisting of a 4-layer film (j = 3 ... 6) sand- 
wiched between two semi-infinite Hubbard models for simple 
cubic lattice in the (001) orientation. Crosses on the vertical 
axes indicate bulk Z values for the left and right substrates. 
Temperature T = 0.02. Lines are drawn as a guide to the eye. 

As the third model system, we study a junction in 
which a finite number of strongly correlated-electron lay- 
ers are sandwiched between two weakly correlated met- 
als. We adopt again a uniform model with respect to 
transfer integrals, i.e., tj = ij = t% = t = 1 for all lay- 
ers. We assign a non-zero but moderate Coulomb energy 
Ul.r = 6 to both metal substrates, whereas in the central 
film we choose a larger Coulomb energy U m . In Fig.f3]we 
show the calculated quasiparticle weight Z for a 4-layer 
film as a function of layer index j in the metallic range 
U m < U c . The calculation was carried out using N = 8 
embedded layers, which comprise the central 4-layer film 
and the two outermost layers of metal substrates on both 
sides. Interestingly, in this thin film Z is very close to 
that of the boundary layers of the semi-infinite metal 
with U — Ur shown in Fig. f2J This rapid convergence 
of Z with increasing film thickness may arise partly from 
the peculiarity of the present model in which the one- 
electron DOS is layer- independent. Thus, there appear 
no finite-size effects such as energy-level discretization in 
the one-electron spectrum at the junction. 

We finally consider a metal/insulator/metal junction. 
Figure Ufa) shows the quasiparticle weight Z of a 4-layer 
film with U m = 16 > U c sandwiched between two metals 
with Ul,r = 6 at temperature T — 0.02. As discussed 
by Helmes et al.^ the metallic states decay within the 
insulating layers, so that Z becomes finite in the film. In 
agreement with these NRG results we find that this pene- 
tration depth within the Mott gap is extremely short. To 
demonstrate this more clearly, we plot in Figs. 4(b) and 
(c) the imaginary part of the layer-dependent Coulomb 
self-energy T,j(ito n ) for the outer and inner layers of the 
film, respectively, as a function of Matsubara frequency. 
For comparison, we also show the imaginary part of the 
impurity self-energy in the bulk simple-cubic crystal with 
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FIG. 4: (color online) (a) Quasiparticle weight Z of a 4-layer 
film with (7 m = 16 sandwiched between two metals with 
Ul,r = 6 as a function of layer index j. Lines are drawn 
as a guide to the eye. (b) and (c) Imaginary part of the im- 
purity self-energy of boundary layers (j = 3, 6) and inner 
layers (j = 4, 5) of the thin film as a function of Matsubara 
frequency for U m = 16 (red solid line). The corresponding 
bulk self-energy at the same Coulomb energy is shown by the 
blue dashed line. 



U = 16. Whereas the bulk self-energy diverges as w n — > 
at this Coulomb energy, the film self-energy tends to 
a finite value because of its contact to the neighboring 
metal layers. At the film surface, the finite value is about 
— 10, indicating bad metallic behavior with a rather short 
electron lifetime and a very small quasi-particle weight 
Z w 0.005. In the second layer, the limiting value of the 
self-energy at low uj n is more than one order of magni- 
tude larger than at the surface, implying correspondingly 
shorter electronic lifetime and lower quasi-particle weight 



Z. Thus, apart from a weak exponential bad-metallic 
tail, the Mott gap of the insulating film is virtually im- 
penetrable. 

Electron transport through correlated-electron sys- 
tems such as oxide heterostructures and molecules is 
emerging as an active field of theoretical studies i 23 ' 24 i 30 
To our knowledge, previous studies considered only non- 
interacting metal leads connected to a central region with 
Coulomb interactions. It would be interesting to extend 
the transport theory to the case of interacting metal leads 
as those shown in Figs. [3] and [U 

IV. SUMMARY 

We have presented an efficient embedding scheme for 
performing DMFT calculations for inhomogeneous lay- 
ered systems such as semi-infinite surfaces and het- 
erostructures. In contrast to previous embedding theories 
based on tight-binding basis functions, the embedding 
potential introduced here is determined from a separate 
DMFT calculation for the bulk substrate materials. It 
therefore incorporates not only the one-electron proper- 
ties but also the many-body effects of the semi-infinite 
systems. The interface region in which local Coulomb 
interactions are treated self-consistently via the layer- 
coupled DMFT also includes the first few layers of the 
actual substrates. As examples, we have presented nu- 
merical results for several surfaces and interfaces of the 
single-band Hubbard model. These results demonstrate 
that the represention of the semi-infinite correlated sub- 
strates in terms of complex energy-dependent embedding 
potentials greatly reduces the numerical effort since only 
a small number of layers needs to be explicitly included 
in the layer-coupled DMFT equation. Thus, the study of 
neutral as well as charged heterostructures involving re- 
alistic strongly correlated multi-band materials becomes 
feasible. 
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